Introduction to the NASA SnowEx Snow Albedo 2023 Dataset
Anton Surunis
2024-04-20
The NASA SnowEx 2023 Snow Albedo Field Campaign Dataset
The NASA SnowEx 2023 Snow Albedo Field Campaign took place in burned and unburned boreal forests around Fairbanks, Alaska. The goal of the campaign was to improve understanding of the spatial, temporal, and process-based variability of snow albedo and the uncertainty of snow albedo measurements across scales in boreal forests. The campaign objectives were to capture snow albedo across scales of snow accumulation and snowmelt with coincident snow albedo from ground-based spectrometer measurements, tower-mounted and drone-based radiation measurements, and airborne AVIRIS-NG overflights across boreal forest disturbance history.
Over five weeks from April 1st to May 5th 2023, several teams visited field sites around Fairbanks and collected spectral measurements over 500m-1km transects capturing snow reflectance and snow albedo over gradients of landscape, topography, and forest disturbance variability. During days with favorable weather/clear sky conditions, teams walked transects in teams of three collecting observations of snow spectra using field spectrometers coincident with hyperspectral aerial and satellite observations from above.
The purpose of this tutorial is to provide an introduction to accessing and using the resulting field dataset. First, a review of background information is provided. Then, we cover how to prepare and access the different data points provided in the dataset. Finally, we provide an example of how to calculate derived statistics from the dataset.
Review of Hyperspectral Data
Incoming solar radiation is either reflected, absorbed, or
transmitted (or a combination of all three) depending on the surface
material. This spectral response allows us to identify varying surface
types (e.g. vegetation, snow, water, etc.) in a remote sensing image.
The spectral resolution, or the wavelength interval, determines the
amount of detail recorded in the spectral response: finer spectral
resolutions have bands with narrow wavelength intervals, while coarser
spectral resolutions have bands with larger wavelength intervals, and
therefore, less detail in the spectral response. 
Surface Reflectance vs Albedo
Hyperspectral data is often captured as either albedo or surface reflectance.
Albedo is the proportion of solar radiation that is reflected by a surface integrated over all incoming solar angles. This is accomplished by taking the ratio of down- and up-facing measurements of hemispherical radiation using a wide (180 degree) lens called a remote cosine receptor (RCR). Albedo is a very important property in calculating land surface energy exchange and snow-mass energy balance.
In contrast, surface reflectance is the proportion of solar radiation reflected over a single or very narrow incoming solar angle (usually 4-8 degrees). Surface reflectance is calculated by taking the ratio of reflected solar radiation from a surface relative and that of a white reference.
White references are usually small panels covered in Spectralon - a highly reflective, near-Lambertian substance that reflects and scatters nearly all incoming light equally in all directions.
Surface reflectance allows us to identify varying surface types (e.g. vegetation, snow, water, etc.) as well as specific qualities of those surfaces (e.g., grain size or grain type in a snowpack). While similar measures, the surface reflectance and albedo of a surface can differ considerably, especially at low solar angles where the angle of direct incident light is far off nadir. Further, since snow reflectance is based on reflected light from a white reference, it is essential that the white reference is kept pristine for accurate measurement of surface reflectance.
Field Spectrometers
Field spectrometers are remote sensing instruments that are carried into the field by operators and used to measure surface reflectance and albedo. Field spectrometers are manufactured by many different companies and come in many more different models using different spectral ranges, attachments, and processing software. While it is difficult to account for these differences without instrument intercomparison studies, it is an important fact to keep in mind when comparing hyperspectral data from different spectrometers.
Data Loading and Description
Let’s load in the data. This is a BIG file so it may take some time…
setwd(dirname(rstudioapi::getSourceEditorContext()$path))
files <- list.files(pattern = ".csv")
refl <- read.csv(files)
as_tibble(refl)## # A tibble: 3,543,003 × 21
## id date instrument site transect type attachment orientation lat
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 2 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 3 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 4 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 5 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 6 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 7 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 8 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 9 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## 10 20230407_… 2023… S2 CARI T2 ssr 8deg down 65.2
## # ℹ 3,542,993 more rows
## # ℹ 12 more variables: long <dbl>, spec_time <chr>, depth <dbl>,
## # depth_alt <dbl>, depth_acc <dbl>, elevation <dbl>, slope <dbl>,
## # aspect <dbl>, tags <chr>, rcr_group <int>, wavelength <int>, value <dbl>The dataset includes field spectrometer measurements of snow reflectance, snow albedo, and up- and down-facing bihemispherical radiance/irradiance measurements along with lots of associated metadata. The column descriptions are as follows:
id: Unique ID number of measurement date: The date of measurement
collection
instrument: Code corresponding to the spectrometer identifier (S1 =
Spectral Evolution; S2 & S7 = ASD FieldSpec4)
site: Code of study site (CARI = Caribou-Poker Creek; DEJU = Delta
Junction, CRMF = Creamer’s Field)
transect: Code corresponding to transect where the measurement was taken
(T1 = burned forest; T2 = forested; T3 = open)
type: The type of spectral measurement as recorded by the note taker
(ssr = snow surface reflectance, albedo = calculated snow surface
albedo, albedo_raw = up and down components of snow surface albedo,
irr_raw = irradiance) attachment: The fiber-optic attachment (8deg = 8
degree optic, 4deg = 4 degree optic, rcr = remote cosine receptor)
orientation: Facing of the fiber-optic attachment (down = down-facing,
up = up-facing)
lat: Latitude of measurement as recorded by the GPS unit
(epsg:4269)
long: Longitude of measurement as recorded by the GPS unit
(epsg:4269)
spec_time: Local date and time of measurement as reported by
spectrometer
depth: Snow depth in cm
depth_alt: Altitude as given by the GPS unit
depth_acc: Accuracy of GPS coordinates as recorded by the GPS unit
slope: Slope of the ground surface in degrees calculated from USGS 3DEP
DEM (10m spatial resolution) using GIS software
aspect: Aspect of the ground surface in degrees calculated from USGS
3DEP DEM (10m spatial resolution) using GIS software tags: Notes taken
by notetaker with discrete notes seperated by “#”
rcr_group: Grouping variable for albedo and irradiance
calculations
wavelength: Wavelength measured by spectrometer
value: Value measured by spectrometer at the given wavelength
Data Preparation
First, we replace -9999 (null) values with NA, set negative values to 0 and convert the date column to the “date” data type.
refl <- refl %>%
mutate(across(where(is.numeric), ~if_else(.x == -9999, NA, .x))) %>%
mutate(value = if_else(value < 0, 0, value)) %>%
mutate(date = ymd(date))Finally, we add some grouping variables to our dataset. These variables are fairly abitrary, but, broadly, transect(s) 1 went through burned forests, transect(s) went through unburned forests, and transect(s) 3 went through open areas. The season variable splits the data into three times spans over the field campaign.
refl <- refl %>%
# Landcover grouping column
mutate(landcover = if_else(site == "CRMF","open",
if_else(transect == "T1", "burn",
if_else(transect == "T2", "forest", "open")))) %>%
# Season grouping column
mutate(season = if_else(date < ymd(20230415), "early",if_else(
date >= ymd(20230415) & date < ymd(20230421), "mid", "late"))) %>%
# Factor season so that it is in the right order
mutate(season = factor(season, levels = c("early", "mid", "late")))
as_tibble(refl)## # A tibble: 3,543,003 × 23
## id date instrument site transect type attachment orientation lat
## <chr> <date> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 2 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 3 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 4 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 5 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 6 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 7 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 8 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 9 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 10 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## # ℹ 3,542,993 more rows
## # ℹ 14 more variables: long <dbl>, spec_time <chr>, depth <dbl>,
## # depth_alt <dbl>, depth_acc <dbl>, elevation <dbl>, slope <dbl>,
## # aspect <dbl>, tags <chr>, rcr_group <int>, wavelength <int>, value <dbl>,
## # landcover <chr>, season <fct>Data Exploration: Measurement Locations
Let’s start exploring the data by mapping our measurement locations.
pts <- refl %>%
filter(type == "ssr" | type == "albedo") %>%
distinct(id, .keep_all = T)
pts.color <- colorFactor(
palette = c("ssr" = "Blue", "albedo" = "Orange"),
domain = pts$type)
leaflet() %>%
addTiles() %>%
addCircleMarkers(
data = pts,
color = ~pts.color(type),
label = ~type,
radius = 1) %>%
addProviderTiles("Esri.WorldImagery")Data Exploration: Snow Surface Reflectance
Let’s take a look at just the reflectance data.
First, we filter our data by snow surface reflectance (ssr) measurements only.
refl_only <- refl %>%
filter(type == "ssr")
as_tibble(refl_only)## # A tibble: 2,280,366 × 23
## id date instrument site transect type attachment orientation lat
## <chr> <date> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 2 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 3 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 4 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 5 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 6 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 7 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 8 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 9 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 10 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## # ℹ 2,280,356 more rows
## # ℹ 14 more variables: long <dbl>, spec_time <chr>, depth <dbl>,
## # depth_alt <dbl>, depth_acc <dbl>, elevation <dbl>, slope <dbl>,
## # aspect <dbl>, tags <chr>, rcr_group <int>, wavelength <int>, value <dbl>,
## # landcover <chr>, season <fct>Some of the measurements have exceptionally high reflectance values, so we filter out measurements where reflectance is too high in the visible range.
refl_only <- refl_only %>%
group_by(id) %>%
filter((!any(wavelength < 750 & value >= 1.2)) %>% replace_na(T))
as_tibble(refl_only)## # A tibble: 2,027,514 × 23
## id date instrument site transect type attachment orientation lat
## <chr> <date> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 2 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 3 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 4 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 5 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 6 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 7 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 8 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 9 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## 10 2023… 2023-04-07 S2 CARI T2 ssr 8deg down 65.2
## # ℹ 2,027,504 more rows
## # ℹ 14 more variables: long <dbl>, spec_time <chr>, depth <dbl>,
## # depth_alt <dbl>, depth_acc <dbl>, elevation <dbl>, slope <dbl>,
## # aspect <dbl>, tags <chr>, rcr_group <int>, wavelength <int>, value <dbl>,
## # landcover <chr>, season <fct>Finally, we plot all the reflectance data in a faceted plot.
wvlim <- c(350,2200)
reflim <- c(1e-5, 1.4)
# Faceted snow surface reflectance plot
plt <- ggplot(refl_only,
aes(x = wavelength, y = value, group = id , color = landcover)) +
# Line plot
geom_line(size = 0.75) +
# Graph colors and axis limits
scale_color_discrete(type = c("darkorange", "darkgreen", "dodgerblue3")) +
scale_x_continuous(limits = wvlim, expand = c(.02, .1)) +
scale_y_continuous(limits = reflim, expand = expansion(mult = c(0, .1))) +
# X and Y labels
xlab("Wavelength (nm)") +
ylab("Snow Surface Reflectance") +
# Formatting
theme_bw() +
theme(
axis.text=element_text(size=8,face = "bold"),
axis.title=element_text(size=8,face = "bold"),
strip.text.x = element_text(
size = 8,
face = "bold"
),
strip.text.y = element_text(
size = 8,
face = "bold"
)
) +
# Facet plot by season and site
facet_grid(season ~ site)
plt
Data Exploration: Snow Surface Albedo
Now, we’ll focus on the albedo measurements specifically.
First, we filter the dataset to include only albedo measurements. Second, we filter out measurements that were deemed low quality by the instrument operators.
albd <- refl %>%
filter(type == "albedo" & (is.na(tags) | !str_detect(tags,"#bad#")))
as_tibble(albd)## # A tibble: 421,596 × 23
## id date instrument site transect type attachment orientation lat
## <chr> <date> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 2 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 3 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 4 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 5 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 6 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 7 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 8 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 9 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 10 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## # ℹ 421,586 more rows
## # ℹ 14 more variables: long <dbl>, spec_time <chr>, depth <dbl>,
## # depth_alt <dbl>, depth_acc <dbl>, elevation <dbl>, slope <dbl>,
## # aspect <dbl>, tags <chr>, rcr_group <int>, wavelength <int>, value <dbl>,
## # landcover <chr>, season <fct>Let’s plot the results (grouped by campaign period, transect, and landcover type) and see how our albedo measurements look.
# Graph limits
wvlim <- c(350,2200)
reflim <- c(1e-5, 1.5)
# Facetted albedo plot
plt <- ggplot(albd, aes(x = wavelength, y = value, group = id , color = landcover)) +
# Line plot
geom_path(linewidth = 0.75) +
# Set colors and axis ranges
scale_color_discrete(type = c("darkorange", "darkgreen", "dodgerblue3")) +
scale_x_continuous(limits = wvlim, expand = c(.02, .1)) +
scale_y_continuous(limits = reflim, expand = expansion(mult = c(0, .1))) +
# Set axis labels
xlab("Wavelength (nm)") +
ylab("Snow Surface Albedo") +
# Formatting
theme_bw() +
theme(
axis.text=element_text(size=8,face = "bold"),
axis.title=element_text(size=8,face = "bold"),
strip.text.x = element_text(
size = 8,
face = "bold"
),
strip.text.y = element_text(
size = 8,
face = "bold"
)
) +
# Facet plot by season and site
facet_grid(season ~ site)
plt
Data Exploration: Snow Surface Albedo Summary
Now, we can try summarizing the results with mean albedo and confidence intervals.
albd_avg <- albd %>%
# Calculate average albedo and upper/lower quartiles
group_by(site, season, landcover,wavelength) %>%
summarise(albedo_avg = mean(value), margin = qt(0.975,df=(length(value)-1))*(sd(value))/sqrt(length(value)), n = length(unique(id))) %>%
ungroup() %>%
# Calculate 95% confidence intervals
mutate(ymax = albedo_avg + margin, ymin = albedo_avg - margin) %>%
# Moving around columns
relocate(c(site, season, landcover), .before = wavelength)Aaaand, let’s check out the results.
# Graph limits
wvlim <- c(350,2200)
reflim <- c(1e-5, 1.15)
# Facet plot of average SSA + 95% confidence intervals
plt <- ggplot(albd_avg, aes(x = wavelength, y = albedo_avg, color = landcover)) +
# Line plot with 95% conifidence "ribbons"
geom_line(linewidth = 0.75) +
geom_ribbon(aes(ymin = ymin, ymax = ymax, fill = landcover),
alpha=0.1,
linetype="dashed",
show.legend = F) +
# Set colors
scale_color_discrete(type = c("darkorange", "darkgreen", "dodgerblue3")) +
# Set axis limits
scale_x_continuous(limits = wvlim, expand = c(.02, .1)) +
scale_y_continuous(limits = reflim, expand = expansion(mult = c(0, .1))) +
# X and Y labels
xlab("Wavelength (nm)") +
ylab("Snow Surface Albedo") +
# Formatting
theme_bw() +
theme(
axis.text=element_text(size=8,face = "bold"),
axis.title=element_text(size=8,face = "bold"),
strip.text.x = element_text(
size = 8,
face = "bold"
),
strip.text.y = element_text(
size = 8,
face = "bold"
)
) +
# Facet plot by season and site
facet_grid(season ~ site)
plt
Broadband Albedo
Broadband albedo (BBA) is the ratio of upward and downward bi-hemispherical reflectance over a specific wavelength range. To calculate BBA, we weight the albedo at each band by the amount of incoming solar radiation (called irradiance) at that band and sum all results over the wavelength range. While albedo is calculated individually over each measurement band, calculations of BBA produce a single value of albedo over a given spectral range. Some common spectral ranges are shortwave BBA (0.25 μm to 5.0 μm), ultraviolet BBA (0.4 μm to 0.7 μm), and visible BBA (0.4 μm to 0.7 μm). Broadband albedo is important for calculating impurities in snowpack, especially ones that absorb light at all short wavelengths such as black carbon.
Data Exploration: Broadband Albedo
Let’s clean up irradiance measurements, apply corrections, and calculate irradiance weights
irr <- refl %>%
# Only use "up"-facing RCR measurements
filter(orientation == "up") %>%
# Remove measurements deemed "bad" by the field team
filter((is.na(tags) | !str_detect(tags,"#bad#"))) %>%
# Filter out noisy NIR and SWIR wavelengths
filter(!(dplyr::between(wavelength, 1350,1420))) %>%
filter(!(dplyr::between(wavelength, 1820,1930))) %>%
filter(wavelength <= 2400) %>%
# Calculate the average irradiance on each day on each transect
group_by(date, instrument, site, transect, wavelength) %>%
summarise(value = mean(value, na.rm = T)) %>%
group_by(date, instrument, site, transect) %>%
# Calculate irradiance weights
summarise(weight = value/sum(value), wavelength) %>%
# Cleaning up
relocate(wavelength, .before = weight)
as_tibble(irr)## # A tibble: 37,380 × 6
## date instrument site transect wavelength weight
## <date> <chr> <chr> <chr> <int> <dbl>
## 1 2023-04-07 S2 CARI T2 350 0.000302
## 2 2023-04-07 S2 CARI T2 351 0.000313
## 3 2023-04-07 S2 CARI T2 352 0.000323
## 4 2023-04-07 S2 CARI T2 353 0.000330
## 5 2023-04-07 S2 CARI T2 354 0.000336
## 6 2023-04-07 S2 CARI T2 355 0.000341
## 7 2023-04-07 S2 CARI T2 356 0.000344
## 8 2023-04-07 S2 CARI T2 357 0.000347
## 9 2023-04-07 S2 CARI T2 358 0.000351
## 10 2023-04-07 S2 CARI T2 359 0.000357
## # ℹ 37,370 more rowsNow we calculate shortwave BBA and clean up the new dataset.
bba <- albd %>%
# Join the irradiance weights to the albedo dataset
left_join(irr, by = join_by(date,instrument,site, transect, wavelength)) %>%
# Apply irradiance weights to albedo measurements
mutate(value = value * weight, .keep = "unused") %>%
# Create a new ID column for later
mutate(id = paste(date, instrument, site, transect, rcr_group, sep = "_")) %>%
# Group and sum across spectrum to calculate BBA
group_by(date, instrument, site, transect, rcr_group) %>%
mutate(bba = sum(value, na.rm = T),
.keep = "unused") %>%
# Drop BBA results that are greater than 1
filter(bba < 1) %>%
# Drop the wavelength column as it is no longer needed
select(-wavelength) %>%
# Remove any duplicate measurements based on the ID we created earlier
distinct(id, .keep_all = T)
as_tibble(bba)## # A tibble: 107 × 22
## id date instrument site transect type attachment orientation lat
## <chr> <date> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 2023… 2023-04-07 S2 CARI T2 albe… rcr <NA> 65.2
## 2 2023… 2023-04-07 S4 CARI T1 albe… rcr <NA> 65.2
## 3 2023… 2023-04-07 S4 CARI T1 albe… rcr <NA> NA
## 4 2023… 2023-04-07 S4 CARI T1 albe… rcr <NA> NA
## 5 2023… 2023-04-07 S4 CARI T1 albe… rcr <NA> NA
## 6 2023… 2023-04-07 S4 CARI T1 albe… rcr <NA> NA
## 7 2023… 2023-04-13 S2 DEJU T1 albe… rcr <NA> 63.9
## 8 2023… 2023-04-13 S2 DEJU T1 albe… rcr <NA> 63.9
## 9 2023… 2023-04-13 S2 DEJU T1 albe… rcr <NA> 63.9
## 10 2023… 2023-04-13 S2 DEJU T1 albe… rcr <NA> 63.9
## # ℹ 97 more rows
## # ℹ 13 more variables: long <dbl>, spec_time <chr>, depth <dbl>,
## # depth_alt <dbl>, depth_acc <dbl>, elevation <dbl>, slope <dbl>,
## # aspect <dbl>, tags <chr>, rcr_group <int>, landcover <chr>, season <fct>,
## # bba <dbl>Let’s make a summary plot.
# Graph limits
bbalim <- c(0, 1)
# Faceted boxplot (BBA)
plt <- ggplot(bba, aes(x = transect, y = bba, fill = landcover)) +
# Boxplot
geom_boxplot() +
# Choose colors
scale_fill_manual(values = c("darkorange", "darkgreen", "dodgerblue3")) +
# Set graph limits
scale_y_continuous(limits = bbalim, expand = expansion(mult = c(0, .1))) +
# X and Y labels
xlab("Transect") +
ylab("Broadband Albedo (350-2400 nm)") +
# Formatting
theme_bw() +
theme(
axis.text=element_text(size=8,face = "bold"),
axis.title=element_text(size=8,face = "bold"),
strip.text.x = element_text(
size = 8,
face = "bold"
),
strip.text.y = element_text(
size = 8,
face = "bold"
)
) +
# Facet plot by season and site
facet_grid(factor(season, levels = c("early", "mid", "late")) ~ site)
plt
Wrapping it up
What we learned:
-Basic intro on how to access and prepare the SnowEx 2023 Snow Albedo
dataset
-How to map points and filter the dataset for specific types of
data
-How to plot and summarize snow surface reflectance and snow
albedo
-How to calculate and summarize broadband albedo